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We study the ground-state properties of a system of identical classical Coulombic point particles, 
evenly distributed between two equivalently charged parallel plates at distance d; the system as a 
whole is electroneutral. It was previously shown that upon increasing d from to cxd, five different 
structures of the bilayer Wigner crystal become energetically favored, starting from a hexagonal 
lattice (phase 1, d — 0) and ending at a staggered hexagonal lattice (phase V, d — >■ oo). In this 

■ paper, we derive new series representations of the ground-state energy for all five bilayer structures. 
I The derivation is based on a sequence of transformations for lattice sums of Coulomb two-particle 

potentials plus the neutralizing background, having their origin in the general theory of Jacobi theta 
r—j • functions. The new series provide convenient starting points for both analytical and numerical 

progress. Its convergence properties are indeed excellent: Truncation at the fourth term determines 
in general the energy correctly up to 17 decimal digits. The accurate series representations are used 
to improve the specification of transition points between the phases and to solve a controversy in 
previous studies. In particular, it is shown both analytically and numerically that the hexagonal 
phase I is stable only at d = 0, and not in a finite interval of small distances between the plates as was 

■ anticipated before. The expansions of the structure energies around second-order transition points 
I ■ can be done analytically, which enables us to show that the critical behavior is of the Ginzburg- 

^ , Landau type, with a mean-field critical index /3 = 1/2 for the growth of the order parameters. 
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PACS numbers: 82.70.-y, 52.27.Lw, 64.70.K-, 73.21.-b 



, I. INTRODUCTION 

. 

O ■ Classical charged particles, confined in a two-dimensional (2D) layer and interacting via the usual three-dimensional 
^ O • Coulomb potential, exhibit a crystallization into a Wigner hexagonal structure, when kinetic energy is small compared 
] to potential energy<ii2 We shall be interested here in bilayer systems, that describe several properties of real physical 
^— I ■ systems in condensed and soft matter, such as semiconductors,'^ quantum dots,'' boron nitrider' laser-cooled trapped 
K*" ] ion plasmas,— dusty plasmas^ and colloids.— For a recent review of numerical methods for quasi-2D systems with 
. long-range interactions, see Refi^ In addition, the creation of a bilayer Wigner crystal on two charged plates at some 
(N- distance is of primary importance in the study of "anomalous" strong-coupling effects such as like-charge attraction 
or overcharging 

In this paper, we study the ground-state properties of a classical one-component plasma of identical Coulombic 
particles of the charge — e, evenly distributed between two plates of the same homogeneous fixed charge density ae 
which are at distance d. The total surface density of the particles is n, the particle density in each layer is ni = n/2. 
The overall electroneutrality of the system is ensured by the condition ni = a. Th e phase diagram of the system at 
temperature T = is determined by a single dimensionless parameter rj = = dy/a. By comparing the static 

energy of various lattices, five distinct phases were detected to be stable (providing global minimum of the energy) in 
different ranges of 77.^^"^^ In order of increasing ij, these phases are: a hexagonal lattice (I) for 77 e [0, rif\, a staggered 
5^ ' rectangular lattice (II) for 77 £ [ij^, 773], a staggered square lattice (III) for 77 £ [ij^, 773], a staggered rhombic lattice (IV) 
, for 77 S [773, 774] and a staggered hexagonal lattice (V) for ij £ [7/4, cxd]; although we use an index c in 77"^, the transition 
point 77"^ from one structure to the other is not necessarily a critical point. The structures are pictured in Fig. [TJ The 
different symbols correspond to particle positions on the opposite surfaces. The primitive translation vectors of the 
Bravais lattice on one of the surfaces are denoted by ai and a2. 

The ground-state structures I, III and V are "rigid" , i.e. they have fixed (77-independent) primary cells within their 
region of stability. The structures II and IV are "soft" , i.e. the shape of their primary cells is varying with increasing 
77, within their region of stability. We now outline the basic characteristics of the structures.—"— 

• Structures I, II and III: Within one single layer, the structure corresponds to a rectangular lattice with the 
aspect ratio A ~ |a2|/|ai|. The equivalent structures on the two layers are shifted with respect to one another 
by a half period, i.e. by (oi + a^)!^. 

Structure I with A = -x/S arises naturally in the simple limit 77 — >■ 0, where the bilayer structure reduces to a 
single layer, which is known to crystallize in a hexagonal (equilateral triangular) lattice i^'^'^^ An open question is 
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Structures I, II, III 



Structure IV 



Structure V 



FIG. 1. Ground-state structures 
I-V of counter-ions on two paral- 
lel equivalently and homogeneously 
charged plates. Open and filled 
symbols correspond to particle posi- 
tions on the opposite surfaces. The 
primitive translation vectors of the 
Bravais lattice on one of the sur- 
faces are denoted by ai and 02. For 
structures I, II and III, we define 
the aspect ratio as A — |ai|/|a2|, 
so that A = \/3 with structure I, 
1 < A < \/3 for structure II and 
A = 1 for structure III. The dashed 
circle for structure IV is a guide to 
the eye, for identifying those points 
that are equidistant to the ion in the 
circle center. For a more detailed 
description of the structures, see the 
text. 



whether phase I (with the fixed aspect ratio A — \/3) exists only at 77 = or is stable also in a finite interval [0, 
with some rj^ > 0. Some numerical calculations indicate very small, but nonzero values of rj^ = 0.006 (Ewald 
techniquel^) and 0.028 (Monte Carlo simulations--). On the other hand, another study for Yukawa bilayers in 
the limit of infinite screening length indicates that rji = 0, so that a buckled phase of type II preempts structure 

I when 1] is small but non vanishingjSl 

Structure II continuously interpolates between the rigid structures I and III. The value of the aspect ratio A 
then changes smoothly from -\/3 at rjf (phase I) to 1 at the transition point 773 to phase III. It is not clear whether 
or not ril, zero or nonzero, is a standard transition point between phases I and II. The transition between phases 

II and III is continuous (of second order). 

• The structure IV is characterized by an angle 9 between primitive cell vectors ai and 02 of the same length a. 
Increasing i], the angle ip changes continuously from 7r/2 at 773 (continuous transition between phases III and IV) 
up to ?74, where it drops to 7r/3. Together with an additional shift between the sub-lattices on the two layers, 
this corresponds to a discontinuous (first order) transition to phase V. 

• The presence of structure V is expected for large enough rj > rjl: At large separation between the layers, the two 
particle sub-lattices are only weakly coupled and so two staggered hexagonal lattices form the stable structure. 

In this paper, we derive new series representations of the energy for all five bilayer structures. The derivation 
is based on a sequence of transformations for lattice sums of Coulomb potentials plus the neutralizing background, 
having their origin in the general theory of Jacobi theta functions. The series has excellent convergence properties, 
which is convenient for numerical investigations, but is also conducive to analytical progress. It will be used to improve 
the specification of transition points between the phases and to solve the aforementioned controversy concerning the 
stability of phase I. In particular, it is shown both analytically and numerically that phase I is realized only at = 0, 
i.e. rii = 0. The expansions of the structure energies around second-order transition points can be done analytically 
which enables us to derive the critical exponents at phase transition points. The critical behavior is of the Ginzburg- 
Landau type,^® with the mean-field critical index j3 — 1/2 for the growth of the order parameters. A preliminary 
account of this work has appeared in Ref.^^ 

The paper is organized as follows. Sec. |ll] is devoted to a detailed derivation of the series representation of the 
energy for structures I-III. The existence of phase I at 77 = only is established analytically and illustrated numerically. 
The second-order phase transition between phases II and III is then described. The energy of structure IV and the 
second-order phase transition between phases III and IV are treated in Sec. IIIIl Sec. IIVI deals with structure V and 
the first-order phase transition between phases IV and V, while our conclusions are finally presented in Sec. |Vl 
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II. PHASES I-III 



Structures I, II and III are treated on equal footings by considering the general case of structure II, see Fig. [T] For 
one single layer, the 2D lattice points are indexed by jai + ka2, where j, k are any two integers (positive, negative or 
zero) and 

ai=a(l,0), 02 = a(0, A) wHh a = (1) 

VctA 

are the primitive translation vectors of the Bravais lattice. The lattice spacing a is determined by the clcctroneutrality 
condition ni = a with the one-layer particle density ni — l/(Aa^). The aspect ratio A is a continuous parameter 
in the interval [l,-\/3]; as was already mentioned, the limiting cases a/S and 1 correspond to the phases I and III, 
respectively. 



A. Energy of phases I-III 



The dielectric constant of the medium is set to unity for simplicity, and we start by a preliminary remark, valid for 
all phases. Our goal is to compute the total electrostatic energy, including particle-particle, particle-plate, and plate- 
plate interactions. The latter two contributions per unit surface can be derived straightforwardly, and respectively 
read Aira^e'^d and —2TTa^e^d. The sum of both, 2TTa^e^d, thus gives 1/2 of the particle-plate energy, and this is why 
in the subsequent analysis, we shall add to the non trivial particle-particle energy one half of the particle-plate energy 
(also referred to as the particle-background term) . The resulting sum provides the full energy of the system. 

The energy per particle E of the bilayer system consists of the intralayer and interlayer contributions, 

E — -Eintra + -Winter- (2) 

We first consider the intralayer contribution. It is well known that lattice sums involving the pair Coulomb interactions 
exhibit infinities which are canceled exactly by the neutralizing background term.— To maintain mathematical rigor, 
we first restrict ourselves to a disk of finite radius R around a reference particle localized at the origin (0,0). The 
interaction energy due to the discrete Wigner crystal is given by 

f Y: ^=r=fxf' f + k'A^<(-y. (3) 

(j,fo)i(0,0) 

Hereinafter, the omission of the lower and upper values for integer indices j, k automatically means a summation from 
—00 to oo. The interaction of the reference particle with the 2D charge background in the disk is expressed as 

2 io |r| 2aA Jo r ^ ' 

-Eintra is the sum of ([3]) plus (|4]). We intend to rewrite i?intra by using the gamma identity 

1 1 dt 



z \^ Jo Vt 



a common procedure in the field. ^'^^ Each term l/(j^ + fc^A^)^/^ in ^ can consequently be written as 

^ ^'=e-*^%-*^^^^ (6) 



y/f + fc2A2 Jo Vi 

As concerns the background contribution (U), the application of the identity ^ to the term 1/r ~ under 
integration leads to 



^/'^ , 1 /■^/'^ . I r dt 2 1 r'^ dt-K 



dr27rr— = / dr2'Kr 



Altogether, we get 



''Jo 



dt 



Jo Vi Jo Vit 



1 - e 



m 



g-t(fl,/a) 



f + k^A^<(^) . (8) 
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Here, the subtraction of unity is due to the absence of the term {j, k) = (0, 0) in the sum ([3]). Having all contributions 
under the same integration, we are allowed to take the limit R/a — ^ cxi, which removes the exponentially small term 
exp[—t{R/a)^] and the disk constraint for lattice indices. Using the definition of the Jacobi theta function with zero 
argumenls^^ 03{q) = J2j 1^' making the substitution tA — > t, we end up with the result 



23/2 



'TT Jo 



dt r 



93(e-*^)e3(e-*/^)-l-T 



We shall repeatedly use the Poisson summation formula 



(9) 



(10) 



The asymptotic behaviors 



2e-* + ■ 



(11) 



follow immediately. We see that the background charge contribution —n/t correctly cancels the i — > singularity of 
the product of two 6*3 functions inside the square bracket in © and the integral converges. 

The Wigner lattices on the opposite layers are shifted with respect to one another by the vector (ai + a2)/2, 
see Fig. [1] To obtain the interlayer contribution to the energy, we first consider the disk of radius R around the 
(perpendicular) image of the reference particle on the opposite layer. The interaction energy of the Wigner crystal is 
given by 



e 

2a 



E 



U-t; +(fc-i)^A2 + (d/a)2 
The interaction with the background charge is described by 

1 e2 



ae 



R 



2aA 



R/a 



drinr- 



+ {d/af 



(12) 



(13) 



Proceeding as in the previous case and taking into account that d/a — "qy A, we find 



23/2, 



dt 



with the Jacobi theta function 6*2(9) = Si ^-^ • It follows from Eq. pO|) that 



so that the integral in Eq. converges, as it should. 
The total energy per particle E reads 



E{A,r)) 




TT' 

1 

t. 



-tA 



)e2{e 



-t/A^ _ ^ 



(14) 



(15) 



(16) 



Note the invariance of i? with respect to the transformation A — > 1/A, which is physically clear from the configuration 
sketched in Fig[T] (label exchange of the two Bravais vectors ai and a2). From a numerical point of view, there are 
two "dangerous" limits: t Q and t — > 00, that jeopardize accuracy. To simplify the integral representation of E, we 
split the range of integration into two parts, [0, tt] and [tt, (X)], and transform the integral over [tt, 00] to the one over 
[0, tt] by using the Poisson formula ([TU]). For the term containing the product of two ^3 functions, the procedure reads 



dt r 

71 



?3(e-*'")03(e-*/'^)-l 



dt 
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dt 

vt 

Tidt' 




1^' 



kfA/t 



n 

1 

t 



r'^ dt 
/o \/i 

Here, going from the third integral to the fourth one, we applied the substitution t' ~ /t. Similarly, we have 



(17) 



dt 

TT yt 



dt 
—i=t 

7T Vt 
°° dt 

" irdt' 



7a 



{tr/' 

dt 



TT 



Vt 



j k 



where the Jacobi theta function 6*4(5) = 7 "^^^ asymptotic behaviors 

1 - 2e"* 



.4(e-)^:^ JV-/(«) + ..., .4(e-). 



(18) 



(19) 



ensure the convergence of the resulting integral. To summarize this paragraph, the total energy (|16p can be rewritten 
as an integral over the finite interval [0, tt] as follows 



£;(A, ri) 



e^V^ 23/2 V^7o ^ 



" dt 



,(e-*^)03(e-*/^)-l-^ 



+e 



j(e-*^)02(e-*/^) - +e-(-'')'/* [04(e-*'^)(?4(e"*/'^) - 1 



(20) 



The exact cancellation of singular terms near t = in this expression for E represents a numerical obstacle that 
should be circumvented. To accomplish the cancellation analytically, we shall consider the series representations of 
Jacobi theta functions and apply to them the Poisson transformation formula (|10p : after subtracting explicitly the 
singular term, the result appears as a series of special functions. In particular, for the first term in the integral (|20)) 
we obtain 



dt 

Vt 



3(e-*'^)^3(e-*/^)-l-^ 



dt 

Vt 



E 



3,k 



-2^ 



dt 



t3/2 



^g-(7rj)V(At)g-(^fc)''A/t _ ^ 



2V^. 



(21) 



The subtraction of the singularity is equivalent to the omission of the term (j, k) = (0, 0) from the summation. Using 
the substitution t' = t/n'^ and introducing the function 



the last expression in Eq. (PT|) can be written as 



r^/'' dt 

/ ^e~^'e-y" for y > 0, 
Jo t 



2E[^3/2(0,jVA) + z3/2(0,j'A)] + 4 ^ z3/2(0,jVA + fc'A)-2V^. 
i=i j,k=i 



(22) 



(23) 
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Finally, performing the above procedure for all terms under integration in (I20[) . we end up with the series representation 



23/2^^ 




OO 



z3/2(0,jVA) + ^3/2(0, j'A)] + 8 ^ z3/2(0,jVA + fc^A) 

j=l j,k=l 



+2^(-l)^- [^3/2(M^JVA) + ^3/2((7^r7f ,j2A)] +4 ^ (-l)-'(-l)'=Z3/2((^r,)2, jVA + fc^A) 
+4^ 23/2(0,772 + (j-l/2)VA + (fc-l/2)2A)-4V^-7rzi/2(0,7?2) I. (24) 

3,k=l J 

The function z,j(x,y) with .t = is related to the so-called Misra functio n^^i^i which was extensively used in single- 
layer lattice summations]^— For our bilayer system with positive 77 we need the more general function (I22p with 
x > 0. The convergence properties of our series can be anticipated from the asymptotic relation Zv{Q, y) ^ e^'^yir"^^ /y 
In numerical calculations, for a given i] we have to find such A* which provides the minimum value of the energy 
p4| . In practice, the series p4)) must be cut at some j, k = M. We document excellent convergence properties of the 
series by considering the single-layer case (A — \/3, '7 = 0) for which the exact^i 

E{V3,0)/{e^V^) = -1.96051578931989165... (25) 

The cut at M = 1, 2, 3, 4 reproduces this exact value up to 2, 5, 10, 17 decimal digits, respectively. A similar accuracy 
is reached in all considered cases. To be extremely accurate, we apply the M = 5 cut everywhere, and use the 
Mathematica software. 

Another advantage of the series representation ([M)) is the possibility of an explicit expansion of the function 
E{l!^* ,r])/{e'^y/n) around the controversial point 77 = and around the critical point 772. As will be shown, this 
requires an analogous Taylor expansion of our z-functions. 



B. Going from phase I to phase II 



We know22i2i that for the single layer, i.e. ?7 = 0, the structure providing the minimum of the energy is the hexagonal 
lattice with A = -y/S- In what follows, we shall investigate the minimum of the energy ([24|) in the neighborhood of the 
point (A — a/S, 77 = 0). We set A = ^/i — e and consider e to be infinitesimally small. To derive a small-e expansion 
of the energy ([M)) . we first perform this task for its series components. From the integral definition it is easy to 
show that the ^-functions under consideration exhibit an analytic (Taylor) expansion in e of the form 



Z3/2(0, j'A) - 23/2(0, j'\/3) + 6^25/2(0, j^Ts) + ieY^7/2(0,j'V3) + ©(e^), 
23/2(0, //A) - ^3/2(0,^7^3) - e^25/2(0,iVV3) ^27/2(0, jVVs) - ^25/2(0, j^/Vs) 



(26) 
(27) 



23/2(0,iVA + fc'A) = z3/2(0,jV%/3 + fc3V3) + e ( - ) 25/2(0, jVVS + fc^Vs) 



^7/2(0,jVV3-|-fc2V3) - ^25/2(0,jVV3 + A:2V3) 



+ 0{e^). (28) 



Similar expansions can be derived for 23/2((7r?7) , j A), 23/2((7r?7) , j /A), etc. Inserting these expansions into 
we obtain 



E{y/?>-e,ri) ^ E{V3,r]) 



+ fiiv)e + f2iv)e'+0{e^), 



(29) 



where the explicit form of the prefactor functions /i(?7) and /2('7) is written in the Appendix. 

Being close to the point (A = \/3, 77 = 0), we are interested in the small- 7; behavior of the functions /i(r/) and f2{v)- 
The corresponding Taylor expansions in powers of 77^ can be performed explicitly, too. The explicit form of _B(\/3, 77) 
is here immaterial. We find that /i(0) = 0, and, with a high precision. 



/i(77) = -0.5833059875 . . .rj'^ + 0(77'*), /2(r/) = 0.0408440789 . . . + 0(77^). 



(30) 



3x10 




FIG. 2. The difference between the dimensionless energies [E{A,rj) — E{y/3,'q)]/{e'^y'n) versus e = ^/3 — A, calculated 
numerically by using (|24|l for two small values of rj: a) r; = 10~^ and b) 77 = 10~^. For the range of aspect ratios chosen, these 
curves are indistinguishable from the analytical prediction — 0.5833 7)^e + 0.0408 e^, stemming from Eqs. (|29p and H30p . The 
values of e* , which provide the energy minimum in the asymptotic limit 77 — )■ according to p2[) . are depicted by the vertical 
dashed lines for comparison. 



For a fixed rj, the extremum of the energy (1291) appears at e* — \/3 — A* given by the stationarity condition 



de 

Namely, 



0^/1(77) +2/2(77)6*. (31) 



V3 - A* = e*(77) - -7^^ = 7.14064 ...rj^+ Oirj^). (32) 

Since d^E{^/3 — ^;'7)|j-g» ^ hiv) > 0; the extremum is the minimum. The result ([32]) tells us that howsoever small 
the dimensionless distance rj is, the buckled structure II with A < a/3 takes place. In other words, the structure 
I exists only strictly at 77^^ = 0. The fact that in the previous works^^^ the structure I was detected also for very 
small positive values of 77 is probably related to extremely small values of the deviation e* c>c 77^ for these t/'s, which 
are "invisible" by standard numerical methods. Like for instance, the structure I border reported in-'^^^ 77^^ = 0.006 
corresponds to e* = 0.00026 . . .. 

In Fig. [21 we present the plots of the difference between the dimensionless energies [i?(A,?7) — E{\/3,r])]/ {e^ ^yn) 
versus e = — A, calculated numerically by using (|24l) for two very small values of 77: a) 77 = 10~^ and b) 77 = 10~^ 
which are well below/above the previous estimate of the phase I threshold,— respectively. Alternatively, using the 
analytical expressions (12^ and ([50)1 leads to the very same data. The nonzero values of e* , which provide the energy 
minima in the asymptotic limit 77 — > according to formula p2p . are depicted by the dashed lines for comparison. 
We see that the energy minima fit well with the expected e* which is clear evidence for the phase I instability. Note 
extremely small values cx 10^^^ — 10^* of the energy difference, in Fig. [21 which justifies the derivation of an accurate 
formula for the Coulombic energy. 

In Fig. [3l the asymptotic relation (|32|) (dashed line) is tested against numerical minimization of the energy (|24|) 
(solid curve) for small and intermediary values of rj, in the logarithmic scale. The agreement is very good, not only 
for small 77, but in the whole range of stability of phase II (it will be shown in the next subsection that phase II border 
is given by 773 — 0.26276 . . .). 



C. Transition between phases II and III 



Going from phase I to phase II is not a phase transition in the usual sense. However, the symmetry of the energy 
E with respect to the transformation A ^ 1/A has the fixed (self dual) point at A = 1 which is the critical point of 
the phase transition from phase II to III. Let us parameterize A as follows A = exp(e). The symmetry A — >■ 1/A is 
now equivalent to e — > — e, i.e. the energy E is an even function of e. The expansion of E around the critical point 




0. 
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0.000 



o.(K)i ' ' 0.01 ' ' 0.1 

FIG. 3. Going from phase I to II: The test of the asymp- 
totic relation ((32} (dashed line) against numerical min- 
imization of the energy (|24p (solid curve) for small and 
intermediary values of r/, in the logarithmic scale. 



, , , , I 

0.0001 0.001 0.01 0.1 

Tl2-ri 

FIG. 4. Transition between phases II and III: The test 
of the asymptotic relation (|37|l (dashed line) against nu- 
merical minimization of the energy (|24|) (solid curve) , in 
the logarithmic scale. 



A = 1 (e = 0) in small deviation e follows from the representation 



+ g2ij]Y+g^{r^Y + 0[e^). (33) 



The explicit form of g2(j]) is written in the Appendix and g^(r\) is not presented due to lack of space, but has been 
derived. The energy (j33p has the Ginzburg-Landau form, e being the order parameter. In contrast to that mean field 
theory, the expression for our energy is exact. 

The critical point is associated with the vanishing of the prefactor to , 

32(77) =0, 77^ = 0.2627602682.... (34) 

The values of 772 obtained in the previous studies were 0.262,^^which is remarkably precise, 0.282& and 0.27i^ The 
functions gi(r\) and 34(77) exhibit the following expansions around the critical 772: 

52(7/) = -0.4620982808 . . . (7/^ - 77) + ©((77^ - 7/)^), 34(77) = 0.1054378203 ... 4- 0{^l - 77). (35) 

The extremum (minimum) of the energy p3p appears at e* ~ A* — 1 given by the condition 



= = 2g2(??)e*+454(??)e*'. (36) 

For 7] < 772 (the "ordered" phase II), we have one trivial solution e* = which however provides the local maximum 
of the energy. There exist two conjugate nontrivial solutions which yield the needed energy minimum; considering 
one of these solutions, we arrive at 

A* - 1 ^ 6* = ^ ^ 1.48031^;^. (37) 

The critical index /3, describing the growth of the order parameter from its zero critical value via e* cx (7^2 — 77)^, has 
the mean field value 1/2. In Fig. |4j in the logarithmic scale, the asymptotic relation (|37|) (dashed line) is compared 
with the numerical minimization of the energy (j24p (solid curve). 

For 77 > 7^2 (the "disordered" phase III), we have the only solution to (|36)) e* = (or equivalently A* = 1) , i.e. the 
rigid phase III is stable, up to a transition to phase IV discussed in the next section. 

The plot of the lattice aspect ratio A* versus 77, obtained by the numerical minimization of the energy ([M)) in the 
whole stability range of the phase II, is pictured by the solid curve in Fig. [5l A* changes from -\/3 at 7/ = to 1 at 
77 = 773. Numerical data of Goldoni and PeetersiS, (open circles) are also presented for comparison. The asymptotic 
relations (|32p for 77 — and (I57|) for 77 — )• 7/2 are also provided, for completeness. 
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FIG. 5. The stability range of phase II: The plot of the lattice aspect ratio A* versus -q, obtained by the numerical minimization 
of the energy H24|) . is pictured by the solid curve. Numerical data of Ref,''* are presented by open circles. The asymptotic 
relations (1321) for 77 — >■ and (|37p for rj ri2 are depicted by dashed curves. 



III. PHASE IV 



In each of the two layers of phase IV (see Fig. [T]), the elementary cell is the rhombus with angle (p between the 
primitive translation vectors 

ai=a(l,0), 02 = a(cos (p, sin 93) with a = , =. (38) 

Vcr sin If 

The lattice spacing a is determined by the electroneutrality condition ni = ct; there is just one particle per rhombus 
of the surface sin ip and so ni = l/(a^ sin(p). The special case of ip — ■7t/2 corresponds to phase III. 



A. Energy of phase IV 



As before, the energy per particle E of the bilayer structure consists of the intralayer and interlayer contributions, 
E = i?intia + ^'intcr- As conccms thc iutralaycr part, the 2D lattice vectors on one layer are indexed with respect to 
a reference particle on the same layer by r(j, k) = jai + ka2, where j, k are any two integers except for (0, 0). The 
square of the lattice vector can be written as 



\v{j, fc)|2 = a2 ^j^2^ 2 jk COS ^) = [{j + kf cos" {p 12) + {3 - kf sin\^/2)\ 



(39) 



This formula represents a kind of "diagonalization" of |r(j, fc)p in indices. If j + fc is an even integer, we introduce 
new indices n = {j + k)/2 and m = (j — fc)/2 covering all integers except for {n,m) ^ (0,0). If j + fc is an odd integer, 
we introduce indices n = (j ' + A; + l)/2 and n = (j — fc + l)/2 covering all integers. Thus the interaction energy due 
to the Wigner crystal can be expressed as 



e 

2a 



E 



1 



j,k 

(j,fc)^(0,0) 



|r(j,fc)| 



e 



E 



1 



„,„"'"(o.o) cos2(^/2)+m2 sin2(<^/2) 



n,m ^{n - 1/2)2 cos2(^/2) + (to - 1/2)2 sm'^{ip/2) 



(40) 



Adding to this expression the interaction with the neutralizing background and using the gamma identity ([5|) in close 
analogy with the previous section, we obtain 



E\„ 



dt 



3(e-*^)03(e-*/'') 



.(e-^)02(e-/^)-^]} 



(41) 
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where 5 = tan{ip/2). 

The Wigner lattices on the opposite layers are shifted with respect to one another by the vector (oi + a2)/2. To 
determine the interlayer contribution to the energy, we first consider the square of the vector between the reference 
particle on one layer and the vertices of the Wigner crystal on the other layer at distance d: 

|r(j, = [{j ^ 1/2)2 ^ _ ^/2)2 + 2{j - l/2)(fc - 1/2) cos^] + 

= a' [U + k-lf cos2((^/2) + (j - kf sm\^/2) + {d/af] . (42) 

Thus the interaction energy with the Wigner crystal reads 



2^^ |r(j, fc)| " 4^ 



E 



1 



n,m \l{n - 1/2)2 COs2(^/2) +TO2sin2(^/2) +d2/(2a)2 

1 



n,m -^7i2cos2((^/2) + (m - 1/2)2 siu^ (^/2) + ^2 /(2a)2 
Adding the background term and using the gamma identity, we find that 



dt 



(43) 



(44) 



The total energy per particle E is the sum of (|4ip and (jH]). Note the invariance of E with respect to the trans- 
formation (5 — ?> 1/S. With respect to the definition oi S = ta,n{(f/2), this symmetry is equivalent to the obvious one 
ip ^ TT — if. Following subsequently similar lines as in previous sections, the integral range [0,c»] can be reduced to 
[0,7r] by using the Poisson formula ([TT 



1 r^i 



-ift/2 



92{e-'')9s{e-'/')~^ 



1-^ 
t 



h{er'')e2{e-*'')-^- 



03(e-")04(e-*/')-l 



e4(e-**)^4(e~*/*) - 1 



(45) 



Applying again the Poisson transformation formula (jlOp to the series representations of Jacobi theta functions, the 
singular f — ^ terms are canceled explicitly and we end up with the representation of the energy per particle E in 
terms of 2;-functions defined in ([2^ : 

= ^1 £[2 + (-l)'] [^3/2(0, jV^) + z„MjH)] + 2 [2 + (-l)^(-l)'=]z3/2(0, jV^ + kH) 

+2 ^ Z3/2(0, (j - 1/2)7-5 + {k- l/2f5) + ^[1 + (-1)^] [z,/,{{7rr^f/2,f/S) + z,/2{{7Tv)V2,fS)] 

oo 

+2 J2 [{-iy + {-ir]zs/2{{7rv)V2,f/S + k'S) 

3.k=l 
oo 

+2 J2 Mo, rf/2 + (j - 1/2)7^ + k'S) + Z3/2(0, vV^ + {j - 1/2)'S + k^S)] 



+ E [^3/2(0' '^VS + U - 1/2)V'5) + 23/2(0, 7,2/2 + (j - 1/2)2,5)] - 3^/^ - ^^1/2(0, vV'^) 



(46) 



B. Transition between phases III and IV 



The symmetry of the energy E with respect to the transformation S ^ 1/d has the fixed point at 5 = 1 which is 
the critical point of the phase transition between the phases III and IV. Parameterizing (5 as (5 = tan((^/2) = exp(— e). 
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the symmetry takes form e — > — e and the energy E is an even function of e. The expansion of E around the critical 
point 5—1 (equivalent to 61 = 7r/2 or e = 0) in small e follows from the representation 

Eie-^ = + + h,ir,)e' + 0{e% (47) 

The explicit form of ^2('7) is presented in the Appendix; The expression for h^^r]) is too lengthy to be given, but is at 
our disposal. 

The critical point is associated with the vanishing of the prefactor of e^, 

h^ir]) =0, ry^ = 0.6214809246.... (48) 

n=V3 

The values of 773 obtained in the previous studies are 0.622^^^ 0.59^ and 0.62.^^ The functions /i2(f]) and /i4(?y) exhibit 
the following expansions around the critical ryg: 

h2{Tj) = -0.2675826391 . . . (r/ - ry^) + 0{{ri - r^lf), hi{ri) = 0.0863245072 . . . + 0(77 - 77^). (49) 

The extremum (minimum) of the energy (j47p appears at e* ~ 7r/2 — Lp* given by the condition 

d E{e-\7^) 



= 2/12(77)6* +4/i4(?7)e (50) 

For rj < Tj^ ("disordered" phase III), we have the only solution e* = (or equivalently ip* — n/2) which provides 
the energy minimum, i.e. the rigid phase III is stable. For ij > rj^ ("ordered" phase IV), the trivial solution e* = 
becomes unstable. The couple of conjugate nontrivial solutions, which provide the energy minimum, implies 



-2-^'-[-im^ == 1.24494^^ (51) 



The critical index /3 has again the mean field value 1/2. In Fig. [Bl in a log-log scale, the asymptotic relation ([5T|) 
(dashed line) is tested against numerical minimization of the energy (|46p (solid curve). In the upper inset, we show 
the dependence of the energy on the logarithm of the angle parameter 6 = tan((/3/2) for 77 = 0.5, where the phase III 
with ^ = 1 is stable. In the lower inset, the analogous plot is presented for 77 = 0.7, where the phase IV with 5 ^ 1 is 
stable; phase III corresponds in fact to a local maximum of the energy. Note the symmetry of the energy with respect 
to the transformation S ^ 1/6 or, equivalently, In 5 — —\nS. 

The plot of the angle parameter d* = tan(95*/2) versus r/, obtained by numerical minimization of the energy (j24p 
in the whole stability range of the phase IV, is displayed by the solid curve in Fig. [T] S* changes from 1 at 77 = 773 
(transition point from phase III to IV) to S'^ = 0.69334 ... at 77 = 774 (transition point from phase IV to V, see the 
next section). Numerical data of Goldoni and Peetersi^ (open circles) are presented for comparison. The asymptotic 
relation ((5T|) for 77 773 is depicted by the dashed curve. 



IV. PHASE V 



In a single layer of the phase V (see Fig. [IJ, the elementary cell of the hexagonal lattice is the rhombus with the 
angle 7r/3 between the primitive translation vectors 

ai=a(l,0), a2 = |(l,\/3) with a = (52) 

The lattice spacing a follows from the electroneutrality condition ni — a; there is just one particle per rhombus of 
surface y/3a'^/2, so that ni = 2/(V3a^). Note that the images of vertices on the opposite layer are localized in the 
center of triangles and not rhombuses, as was the case of phase IV. There is no continuous way to pass from phase 
IV to phase V. 




Tl-ri, 



FIG. 6. Transition between phases III and IV: The test 
of the asymptotic relation (|51[l (dashed line) against a 
numerical minimization of the energy ((46]) (solid curve), 
in the logarithmic scale. The content of upper and lower 
insets is commented in the text. 




0.75 



FIG. 7. Stability range of phase IV: The plot of the 
angle parameter 5* versus -q, obtained by the numerical 
minimization of the energy H46p . is shown by the solid 
curve. Numerical data of Ref.>lS are presented by open 
circles. The asymptotic relation (|5ip for rj ^ rj^ is de- 
picted by the dashed curve. 



A. Energy of phase V 



To study the intralayer contribution to the energy of the reference particle at the origin, we first consider the 
Wigner crystal of lattice vectors r(j, k) = jai + ka2, where integers (j, k) ^ (0, 0). The square of the lattice vector is 
expressible as 



|r(j,fc)|^ 



{f + k'+jk)^^[3{j + k)' + {j^kr 



(53) 



In analogy with phase IV, we introduce new n, m indices for each of the cases j + k being an even and odd integer. 
The interaction energy due to the Wigner crystal then reads 



e 

2a 



E 



1 



|r(j,fc)| 



e 

2a 



E 



1 



E 



^3(71 - 1/2)2 + (m- 1/2)2 



(54) 



.(,1,771)^(0,0) 

Adding to this expression the interaction with the neutralizing background and using the gamma identity, we find 



intra 



1 



dt 



{ [03(e-*/^)^3(e-*^) - 1 - + [e2(e-*/^)02(e-*^) -j]}. 



(55) 



The hexagonal lattices on the opposite layers are shifted with respect to one another by the vector (ai + a2)/3; 
note that the factor 1/3 differs from 1/2 of the previous phases I-IV. To determine the interlayer contribution to the 
energy, we first consider the square of the vector between the reference particle on one layer and the vertices of the 
Wigner crystal on the other layer at distance d: 



|r(j, fc)|2 = a2 ^ ^/3)2 ^ ^ ^/3)2 ^ ^ 1/3)^^ ^ ^^3)] + d" ^ - [3(j + k + 2/3)^ + [j ~ kf] + d\ (56) 

Going from (j, fc) to integers (n,m), the interaction energy with the Wigner crystal on the opposite side takes the 
form 



E' 



1 



2aj^\vij,k)\ 



e 

2a 



E- 



1 



E' 



1 



v/3(n + 1/3)2 + to2 + (rf/„)2 ^3(7^ - 1/6)2 + (m - 1/2)2 + (rf/„)2 

Adding the background term, using the gamma identity and the readily derivable relations 



E' 



-3*0 + 1/3)" 



1 r 

2 



(e-*/3)_^3(,-3*) 



E' 



,-3t(i-l/6)" 



(57) 



(58) 
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we find that 



dt / 1 



73 



{[03(e-*/^)^3(e-*^) 



(59) 



The total energy per particle E is given by the sum of ([55]) and ([5^ . The Poisson formula enables us to reduce 
the integral range to [0, tt], 



E{r]) 1 r dt 



1 .,2 



V3 



-3ift/2 



+ 1 



V3 



-3rft/2 



2(e-*/^)e2(e-*^) - J 



+ I 1 _ le-^-'')'/^^*) + ^^g-3MV(2t) j [04(e-*/v/3)0,(e-*^) - 1 



(60) 



In terms of the functions 

r dt 
Jo Vt 



e2(e-*/^)02(e-*^) - ^ 



2^(-l)J" [z3/2(a;,y + jVV3) + Z3/2(x,y+/V3)] +4 ^ (-l)^(-l)'=Z3/2(a:, y + //^^ + ^'VS), (61) 



J3(x,y)EE^^-^e--*/-%-^-^/* 



'3(e-*/^)03(e-*^)-l-y] 

OO 

2^ [z3/2(a:,y + jV%/3) + 23/2(a:,y+j2y3)] + 4 ^ zy2{x,y + f / + ^Vi) - ^z^/^{x,y), (62) 
—1 i,fc=i 



^4(2;, y) 



r dt 
Jo Vt 



4(e-*/^)04(e-*^) - 1 



= 4 ^ zy2{x, y + {j~ l/2)V\/3 + (fc - 1/2)2^3) - 7rzi/2(x, y), 



i,fe=i 



the energy per particle E is expressible as 



+ 



2/3(0, 0) - i/3((^^)V2, 0) - ^/3(0, 77V2) + ^/3(3(^r,)V2, 0) + ^13(0, 3r;V2) 



/2(0,0) - i/2((7r,7)V2,0) + ^/2(3(7rr;)V2,0) 



/4(0, 0) - i/4(0, r;V2) + ^14(0, 37?V2) 



(63) 



.(64) 



B. Transition between phases IV and V 

Increasing r] from 773, phase IV is stable up to the point rjl at which the energy of phase IV (|46p . evaluated at 5* 
which minimizes this energy, equals to the energy of phase V (|64p . Our result is 



ril = 0.73242.... 

The values of 77^ obtained in the previous studies were relatively dispersed: 0.732)^^ 0.70^ and 0.87i^ 



(65) 
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The phase transition is of first order since the energies of phases IV and V have as functions of rj different slopes 
which causes the discontinuity of the first derivative of the energy with respect to rj at the transition point rj^. The 
angle parameter i5, which minimizes the energy of the phase IV at the critical point r/l, is found to be 6'^ = 0.69334 . . .. 
Since 6 = tan((p/2), the corresponding angle — 69.4702 . . .°; this angle is very close to the estimate ip'^ ~ 69.48° of 
the workji^ Going from phase IV to V, the angle skips to 60° as is indicated in Fig. [T] 



C. Discussion 



Two different "sum-rules" can be derived, that allow for a critical assessment of the results obtained. The simplest 
one relies on the geometrical proximity between structures I (a single hexagonal crystal) and V (two hexagonal crystals 
at half density). For large distances, the two crystals decouple and we have, making use of straightforward notations, 

Ei{V3,r] = 0) = V2Ev{V'^ 00). (66) 

With our series representations for Ej and (IM)) for Ey, this identity holds. Another more subtle constraint 
follows from a combination of elementary geometric considerations)^^ which impose that 

Eviv^O) ^^^^Ej[V3,v^0)- (67) 

Wc have also checked that this identity holds with the expressions provided above (note though that 77 = lies outside 
the stability range of structure V) . 

Finally, it is interesting to consider both the large small and large distance behavior of the energy. For small rj, it 
can be shown that both structures I and II share the same energy expansion, up to order 77'^ included: 

Eii{A*,rj)^Ei{Tj)+0{7j^) (68) 

where A* is the previously introduced optimal aspect ratio that minimizes Ejj{A* , 77) for a given 77. Explicit calculation 
up to order 77^ shows that 

^ + _ 2.59372 . . . + 0(,3), (69) 

where the precise value of Eii{V3,0) = Ei{0) has been given in Eq. (j25p . We note that the linear term in (j69p 
generates a contact pressure — 27rcr^e^ for 77 — ;> 0. A similar term was reported where however the term in r]^ 
differs from ours by a large factor (0.2122 instead of 2.5937) 

At large distances, the relevant phase is structure V, from which the inter-platc pressure follows. The large 77 case 
is encoded in the small t limit of Eq. (1501) , or equivalently, Eq. ([5^ . A saddle point argument leads to 

We recover an expression already obtained iuji^ at variance with other approaches.™ Taking the 77 derivative and 
remembering that n ~1a yields the inter-plate pressure 

The ?7 dependence is well known, since it can be written exp(— God), with Gq the modulus of the first reciprocal lattice 
vector. It should be noted though that the prefactor differs from the often reported 27r(cre)^ (see e.gJ^), by a factor 
3. 



V. CONCLUSION 



The system of classical charged particles, forming a sequence of bilayer Wigner structures at zero temperature as 
the distance between the plates is increasing, has a rather long history. We have presented here a new method to 
calculate the Coulomb ground-state energy of each Wigner structure. Based on a series of transformations and using 
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FIG. 8. Summary of phase transition scenario. The rounded off values of the thresholds are mentioned. Structure I is realized 
at 77 = only (vanishing inter-plate separation). 

general properties of the Jacobi theta functions, we expressed the energies in terms of quickly converging series of 
the functions z^{x,y) defined in (I22p . The presence of the neutralizing background manifests itself simply as the 
subtraction of singularities of the Jacobi theta functions under an auxiliary integration. 

Numerical evaluation of the series requires modest computer and programming facilities, and at the same time 
provides extremely accurate estimates of the energy. We took advantage of this feature, supplemented by analytical 
work, to improve and complete previous studies in three aspects: 

• There was a relatively large dispersion in the determination of the transition points between phases; this concerns 
especially the first-order phase transition between phases IV and V. Our method improves significantly the 
location of all transition points, that can be worked out with arbitrary precision. Figure [S] gives an overview of 
the sequence of phases together with the corresponding thresholds. 

• We resolved, analytically and numerically, a previous controversy about the stability phase I, thereby corrob- 
orating the findings of Ref.-^ We found that this phase is stable only at zero distance between the plates, 
rj = Q = d. To confirm numerically this result, we worked with extremely small values of the energy differences 
oc 10"""^^ — 10~® for distances 77 = 10"'^ and 10~^ (see Fig. [2]), which are "invisible" by standard numerical 
methods. The agreement between the 77 — asymptotic relation calculated analytically by using the 
Taylor expansions of the functions z^{x, y), and the numerical minimization of the energy, presented in Fig. |3l 
is excellent. 

• The expansions of the structure energies around second-order transition points can be done analytically which 
enables us to specify the critical phenomena at the phase transition points; see the expansions pertaining to 
the transitions from phase II to III p7p and from phase III to IV ([5T|) . The agreement between these analytic 
formulas and numerical minimization of the ground-state energy is very good; It can be appreciated in Figs. 0] 
and[Sl Quite expectedly for a zero temperature situation, the critical behavior is always of the Ginzburg-Landau 
type, with the mean-field critical index /3 = 1/2 for the growth of the order parameters in the "ordered" phases. 

It is clear that our method can be directly generalized to other problems concerning the lattice summations over 
pair interactions, not only the Coulomb ones. The bilayers with repulsive Yukawa interactions, extensively studied in 
the past^i2i or with inverse power laws^^ deserve our attention. We additionally emphasize that the ground states 
under consideration here are such that the ions are distributed evenly (50% on each plate): in other words, the ionic 
surface density of one Wigner crystal on a given plate is a, and coincides with the plate homogeneous surface density. 
When dealing with asymmetric plates, this local neutrality assumption should be relaxedf^i^ still enforcing global 
neutrality. This makes the asymmetric problem significantly more complex, and an interesting perspective for future 
work. Finally, consideration of dielectric jumps between the walls and the interstitial slab is also a relevant venue for 
forthcoming investigations. 
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APPENDIX 



The prefactor functions /i (77) and /2 (77) of the expansion (1^^ read 




+4E 



1 



25/2(0, + (j - 1/2) VV3 + ik- 1/2)2 V3) I 



/2(r?) 



23/2 

00 



4E T^7/2(0'J'^) + l^Zj/2{0,f/V3) ~ -|^Z5/2(0,j2/V3) 



L 

33/2 



,.2 r 



27/2(0,j2/V3 + fc2V3) - J^Z5/2iO,f/V3 + k^VS) 



y27/2((7r77)2,/y3) + ^zy.Hnrjf , f / VS) - ^z^Minvf ,f /VS) 



+4 E (-iF(-i)' 



+4E 



1\^ 1 



27/2(0, + U - l/2)2/V3 + (fc - 1/2)2^3) 



"3^ (j-^) 25/2(0, 7?2 + (j-l/2)2/V3+(fc- 1/2)2^3) 
The prefactor function 52 {v) of expansion ([33|) takes the form 



^|2E[A7/2(0,/)-A5/2(0,J^)] 



+ 2 E [(j'-fc')'2,/2(0,J^ + fc^)-(/+fc^)25/2(0,j2+fc2)] 



+ E(-1)' [j'^7M{nrjf,f)-fz,/,{{nrj)\f)] 



+ E (-l)^(-l)'= [if - krz7M{nv)\f + k') - if + k')z,;,{{nrj)\f + k')] 



00 



+ E ( [U 1/2)^ - (fc - 1/2)2] ^^^^(0^ ^2 ^ _ ^/2)2 + (fc - 1/2)2) 
- [ij 1/2)2 ^ _ 1/2)2] ^^^^(0^ ^2 ^ _ 1/2)2 ^ _ 1/2)2)) |_ 
Finahy, the prefactor function /i2('7) of the expansion (H71) can be written 
/^2(r7) 



2rl^£t-^'''^/2(0,/)-J^25/2(0,J^)]+E(-l)qA7/2(0,j2)-j2.5/^ 

00 

+ 2 E [(j'-fc'f2,/2(0,j2 + fc2)-(/+fc2),5/,(0,,2+A;2)] 

+ E (^l)^(-l)'= [if - k'f ^7/2(0, f + k') if + k')z,/,{0,f + fc2)] 
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+ J2 ( [ij - 1/2)' -{k- l/2f]zy2{0, U - 1/2)' + {k- 1/2)2) 

3,k=l 

- [{j - 1/2)2 ^ _ 1/2)2] ^^^^^(0^ _ 1/2)2 + _ 1/2)2)^ 

OO 

+2 5] (-1)^- [(f _ fc2)2^^^^((^^)2/2^_^-2 ^ ^2) _ ^ k'),^^^{{7r^)y2,f + k')] 

j,fe=l 

oo 

+ ^ [(i - 1/2) V/2(0, 77V2 + (i - 1/2)') - (j - 1/2)2^5/2(0, 772/2 + (j - 1/2)2)] 

J=l 

00 

+2 ^ ( [{j - 1/2)2 _ ^2] 2 ^^^^(0^ ^2/2 + _ 1/2)2 ^ ^2) 

j,k=l 

- [{j - 1/2)2 ^ ^2] ^^^^(0^ ^2/2 + _ 1/2)2 ^ ^2)\ I ^ (A4) 
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